function obj = cnt(varargin) 
%
%[StrOutput,info] = MakeCnt(OPTIONS) computes the coordinates of each atoms of a Carbon
%Nanotube with a generic chiral vector and lenght and write it in a xyz
%file.
%m,n = chiral indexs
%ao  = lattice costant
%L   = CNT lenght. The units is specified by the options.
%OPTIONS consists of a pair values and can be omitted.
%MatOutput conteins the coordinates of the atoms in "MatStructure" format. 
%OPTIONS:
%       n         = <5>                         : n chiral index
%       m         = <5>                         : m chiral index
%       a0        = <0.77>                      : Lattice constant
%       L         =  1                          : Lenght in units speified
%       by Units
%       Units     = <PL        |<Angstrom>      : PL stands for principal layer and means
%                                                 the shortest periodic lenght
%       Plane     = <false>    | true           : If true the CNT preparation is showed
%       Preview   = <false>    | true           : If true a fast visualization of CNT is done
%       Filename  = <cnt_n_m>  |                : Name of the output xyz file
%       Sorting   = <true>     | false          : If true the sorting along
%       the growth direction is done
%
%       Position  = <0,0,0>    | x0,y0,z0       : Coordinates of the origin of the CNT
%       Direction = <0,0,1>    | dirX,dirY,dirZ : Direction of CNT growth
%       Capping   = <false>    | true           : Do capping if a cap model is availible 
%       Verbose   = <false>    | true           : If true some informations
%                                                 about the CNT are displayed 
%NOTE: the option pair values can be inserted in any order.
%OUTPUT:
%StrOutput = CellArray containing atoms coordinates and their species.
%info = information about the cnt: info.TotalLenght (in Amstrong) 
%                                  info.diameter (in Amstrong)
%Example:
% MakeCnt('n',12,'m',10,'L',2)
%
% Author: Giuseppe Romano.
% Email: giuseppe.romano@uniroma2.it



%% Parse input
p = inputParser;
p.addRequired('n',@(x)x>=0);
p.addRequired('m',@(x)x>=0);
p.addRequired('L',@(x)x>=0);

%p.addParamValue('help',false,@islogical);
p.addOptional('a0',1.42,@(x)x>0);
%p.addParamValue('L',1,@(x)x>=0);
%p.addParamValue('n',5,@(x)x>=1 && mod(x,1)==0);
%p.addParamValue('m',5,@(x)x>=1 && mod(x,1)==0)
p.addOptional('Units','Angstrom',@(x)any(strcmpi(x,{'Angstrom','PL'})))
p.addOptional('Preview',false,@islogical)
p.addOptional('Plane',false,@islogical)
p.addOptional('Sorting',true,@islogical)
p.addOptional('Position',[0,0,0],@(x) isequal(size(x),[1,3]))
p.addOptional('Direction',[0,0,1],@(x) isequal(size(x),[1,3]))
p.addOptional('Verbose',false,@islogical)
p.addOptional('Capping',false,@islogical)
p.addOptional('FileName','a256t',@ischar)
p.parse(varargin{:});

do_help = 0; %p.Results.help;
m = p.Results.m;
n = p.Results.n;
L = p.Results.L;
a0 = p.Results.a0;
Units = p.Results.Units;
Plane = p.Results.Plane;
Verbose = p.Results.Verbose;
Preview = p.Results.Preview;
pos = p.Results.Position;
dir = p.Results.Direction;

FileName = p.Results.FileName;
if (FileName == 'a256t')
   FileName = ['cnt' '_' num2str(n) '_' num2str(m) '.xyz'];
end


%%Further checks

if (m>n)
   error('Condition 0<|m|<n failed')
end
%% Hexagon parameters
beta = 2 * a0 * 0.8660254; %aphotema times 2
alpha = 3 * a0;

%% Calculate the layer lenght

if (n==m)
   G=-1;
else
   H=(m+n)/(m-n);
   G=(1-3*H)/(1+3*H);
end
    
ss=0; 
%is_int=0;
max_it=1000;
nit=0;
rest = 1;
%while (is_int==0 && nit<max_it)
while (rest > 1e-6 && nit<max_it)
     nit = nit+1;
     ss=ss-1;
     kk=G*ss;  
     rest = abs(round(kk)-kk);
     
     %if (abs(round(kk)-kk)<1e-4)
     %    is_int=1; 
     %end
end

if (nit==100)
   disp('No periodicity condition reached') 
   return       
end

Xpl(1)=(kk+ss)*(alpha/2);
Xpl(2)=(kk-ss)*(beta/2);

Lpl = norm(Xpl);



%% Compute the chiral vector and its ortogonal vector

%Chiral vector%%%%%%%%%%
chi(1)=alpha/2*(m+n);
chi(2)=beta/2 *(m-n);
modChi=sqrt(chi(1)^2+chi(2)^2);
phiChi=atan(chi(2)/chi(1));

%xmax = max(m,n);
%xmin = min(m,n);
chiral_angle = atan((sqrt(3) * m)/(2*n + m)) * 180/pi;
%Length vector
vec=cross([chi,0],[0,0,-100]);
Lvec=vec(1:2);
Lvec=Lvec./(sqrt(Lvec(1)^2+Lvec(2)^2))*Lpl;
%%%%%%%%%%%%%%%%%%%%%%%

%% Find the outer parallelogram where the check must be done

ma=beta/alpha;
Q1(1)=0;
Q1(2)=0;
ma1=-ma;
Q2(1)=Lvec(1);
Q2(2)=Lvec(2);
ma2=ma;
Q3(1)=Lvec(1)+chi(1);
Q3(2)=Lvec(2)+chi(2);
ma3=-ma;
Q4(1)=chi(1);
Q4(2)=chi(2);
ma4=ma;

%Intersection points
P1(1)=(ma1*Q1(1)-Q1(2)-ma2*Q2(1)+Q2(2))/(ma1-ma2);
P1(2)=ma1*P1(1)-ma1*Q1(1)+Q1(2);

P2(1)=(ma2*Q2(1)-Q2(2)-ma3*Q3(1)+Q3(2))/(ma2-ma3);
P2(2)=ma2*P2(1)-ma2*Q2(1)+Q2(2);

P3(1)=(ma3*Q3(1)-Q3(2)-ma4*Q4(1)+Q4(2))/(ma3-ma4);
P3(2)=ma3*P3(1)-ma3*Q3(1)+Q3(2);

P4(1)=(ma4*Q4(1)-Q4(2)-ma1*Q1(1)+Q1(2))/(ma4-ma1);
P4(2)=ma4*P4(1)-ma4*Q4(1)+Q4(2);

%% Compute the largest honeycomb vectors that belong to the outer parallelogram

modV=sqrt((alpha/2)^2+(beta/2)^2);
m_min=0;
m_max=ceil(norm(P2-P1)/modV);
n_min=-ceil(norm(P1)/modV);
n_max=ceil(norm(P4)/modV);


%% Plot chirality

if (Plane)

    figure
    hold on
    axis equal
    title=strcat('CHIRALITY (',num2str(m),',',num2str(n),')');
    set(gcf,'Name',title) 
   
    plot([0,chi(1),chi(1)+Lvec(1),Lvec(1),0],[0,chi(2),chi(2)+Lvec(2),Lvec(2),0],'r') %original inner parallelogram
    plot([0,alpha/2],[0,beta/2],'b','LineWidth',2)  %a1
    plot([0,alpha/2],[0,-beta/2],'b','LineWidth',2) %a2
    plot([P1(1),P2(1),P3(1),P4(1),P1(1)],[P1(2),P2(2),P3(2),P4(2),P1(2)],'g') %original outer parallelogram
    
    for km=m_min:m_max
      for kn=n_min:n_max
          x(1) = alpha/2 *(km+kn);
          x(2) = beta/2  *(km-kn);
          plot(x(1),x(2),'Marker','o','LineStyle','none','Color','b') 
          x(1) = x(1) + a0;
          plot(x(1),x(2),'Marker','o','LineStyle','none','Color','b')
      end
    end

end

%% Construct the rotate inner bound box
bound = zeros(2,2);
bound(1,1)=-modChi/2;
bound(1,2)=modChi/2;
bound(2,1)=0;
bound(2,2)=Lpl;

%% Build the translate vector to origin
x0(1,1)=-modChi/2;
x0(2,1)=0;


%% Build the rotation matrix to z-positive direction
teta=-phiChi;
rot=zeros(2,2);
rot(1,1)=cos(teta);
rot(2,1)=sin(teta);
rot(1,2)=-rot(2,1);
rot(2,2)=rot(1,1);


%% 1PL CNT bilding

C=modChi; %Circumference
D=C/pi;   % Diameter
epsb = 1e-4;

k=0;
for km=m_min:m_max
      for kn=n_min:n_max
          
          %First base point
          
          x(1) = alpha/2 *(km+kn);
          x(2) = beta/2  *(km-kn);
          xrot = rot * x' + x0;
          if  (xrot(1)>=bound(1,1)-epsb && xrot(1)<=bound(1,2)-epsb )
              if (xrot(2)>=bound(2,1)-epsb  && xrot(2)<=bound(2,2)-epsb )
                                  
                 xrot(3)=D/2*(1-cos(2*xrot(1)/D));
                 xrot(1)=D/2*sin(2*xrot(1)/D);                  
                 
                 k=k+1;
                 prin_layer(k,1)=xrot(1);
                 prin_layer(k,2)=xrot(2);
                 prin_layer(k,3)=xrot(3);
                
              end
          end
          
          %Second base point
           
          x(1) = alpha/2 *(km+kn) + a0;
          x(2) = beta/2  *(km-kn);
          xrot = rot * x' + x0;
          if  (xrot(1)>=bound(1,1)-epsb  && xrot(1)<=bound(1,2)-epsb )
              if (xrot(2)>=bound(2,1)-epsb  && xrot(2)<=bound(2,2)-epsb )
                 
                 xrot(3)=D/2*(1-cos(2*xrot(1)/D));
                 xrot(1)=D/2*sin(2*xrot(1)/D);                  
                 
                 k=k+1;
                 prin_layer(k,1)=xrot(1);
                 prin_layer(k,2)=xrot(2);
                 prin_layer(k,3)=xrot(3);
                             
              end
          end
          
      end
end
          

%% CNT building 

if (strcmpi(Units,'PL'))
    Npl=L;
    L_res=0;
else
    Npl=fix(L/Lpl);
    L_res=L-Lpl*Npl;
end

%Replicate along y-positive direction
natom_PL=k;
temp = [];
for kpl=1:Npl
    for k_loc=1:natom_PL
        temp((kpl-1)*natom_PL+k_loc,:)= prin_layer(k_loc,:) + [0,(kpl-1)*Lpl,0]; %
    end
end
 
%Add the last part
for k_loc=1:natom_PL
     if (prin_layer(k_loc,2)<L_res)
         temp(Npl*natom_PL+k_loc,:)= prin_layer(k_loc,:) + [0,Npl*Lpl,0]; 
     end
end

  
for ka=1:size(temp,1)
    %Center to the origin
    temp(ka,3)=temp(ka,3)-D/2;
end

CntAtoms = size(temp,1);
CntLenght = Lpl*Npl+L_res; %Angstrom

%Sorting along y direction

if p.Results.Sorting
   [Y I] = sort(temp(:,2));
   temp2 = zeros(CntAtoms,3);
   for k=1:CntAtoms
      temp2(k,:) = temp(I(k),:);
   end
   temp = temp2;
end


%% Do capping
CappingAtoms = 0;
if p.Results.Capping

    %Search the database for available caps
    DBCAP=importdata('commons/CntCapDB.dat');
    for k=1:size(DBCAP.textdata,1)
        
        if ((str2num(DBCAP.textdata{k,1}) == m)...
            & (str2num(DBCAP.textdata{k,2}) == n)) 
            info.capping = true;  
            namefile = DBCAP.textdata{k,3};
            theta = DBCAP.data(k,1);
            phi = DBCAP.data(k,2);
            gamma = DBCAP.data(k,3);
            
        end
    end
       
    if info.capping
       ful = NBP.mat_structure;
       completename =['commons/FullereneLibrary/' char(strtok(namefile)) '.cc1'];
       NBP.interface(ful,'FileReader',completename);
       sel = NBP.interface(ful,'Selection','All');
      
       theta2 = 90 - theta; 
       phi2   = -phi;
       %Rotation
       NBP.interface(ful,'Transformation','rotation',sel,zeros(1,3),[0 0 1],-theta2); 
       NBP.interface(ful,'Transformation','rotation',sel,zeros(1,3),[1 0 0],-phi2); 
       NBP.interface(ful,'Transformation','rotation',sel,zeros(1,3),[0 1 0],gamma);
       %Splitting
       
       Di = 0;
       middlep = NBP.interface(ful,'Selection','Threshold','y','greater',1e-5);
       NBP.interface(ful,'Transformation','translate',middlep,[0  Di+CntLenght-min([a0/2,CntLenght]) 0]);
       middlem = NBP.interface(ful,'Selection','Threshold','y','less',1e-5);
       NBP.interface(ful,'Transformation','translate',middlem,[0 -min([a0/2,CntLenght])-Di 0]);
       %%%%%%
     
       %NBP.interface(fullerene,'PlotAtoms','MolViewer');
       sel = NBP.interface(ful,'Selection','All');
       CellStructure = NBP.utility(ful,'GetAtoms',sel);
       for ka = 1:size(CellStructure,1)
           temp = [temp;CellStructure{ka,1}];
       end
       CappingAtoms = ful.natoms;
    end
   
end





%% CNT rotate and translate along the user given direction and position
%Rotate

dir=dir./norm(dir);
vec=[0,1,0];
u=cross(dir,vec);
teta=acos(dot(vec,dir)/norm(vec));
    
cosa = cos(teta);
sina = sin(teta);
vera = 1 - cosa;
x = u(1);
y = u(2);
z = u(3);
rot3D = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; ...
          x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; ...
          x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';

  
for ka=1:size(temp,1)
    %Roto-translate
    temp(ka,:)=rot3D*temp(ka,:)'+pos';
end


%% Check Preview option

if Preview

   figure
   hold on
   
   title=strcat('CNT (',num2str(m),',',num2str(n),')');
   set(gcf,'Name',title) 
   
   for na=1:size(temp,1)
       plot3(temp(na,1),temp(na,2),temp(na,3),'Marker','o','LineStyle','none','Color','r'); 
   end

   xlabel('X')
   ylabel('Y')
   zlabel('Z')
   axis equal
   
   dir=dir.*10;
   %plot3([0,dir(1)],[0,dir(2)],[0,dir(3)])
   
   ROTATE3D on;
end



%% Write atom data
sp = 'C';
for ka=1:size(temp,1)
    atpos(1)=temp(ka,1);
    atpos(2)=temp(ka,2);
    atpos(3)=temp(ka,3);
    StrOutput{ka,1}=atpos;
    StrOutput{ka,2}=sp;
end



%% Fill info

info.lenght = CntLenght;
info.diameter = D;
info.start_point = pos;
info.orientation = dir;
info.capping_atoms = CappingAtoms;
info.cnt_atoms = CntAtoms;
info.natoms = CntAtoms + CappingAtoms;
info.chiral_angle = chiral_angle;

%% Write output on the screen

if Verbose
    disp('CNT description')
    disp(' ')
    disp('Total lenght')
    disp(Lpl*Npl+L_res)
    disp('PL Length')
    disp(Lpl)
    disp('Number of PL ')
    disp(Npl)
    disp('Residual Lenght')
    disp(L_res)
    disp('Number of atom per PL')
    disp(natom_PL)
    disp('Chiral Angle')
    disp(info.chiral_angle)
end  


%Building the structure

natoms = size(StrOutput,1);
obj = NanoObj;
for ka=1:natoms
    obj.add_atom(StrOutput{ka,1},'C');
end



end


%fid = fopen(FileName,'w');
%natoms=size(StrOutput,1);
%count = fprintf(fid,'%5g\n',natoms);
%count = fprintf(fid,'%5s\n',' ');
%for ka=1:natoms
    %count = fprintf(fid,'%5s',char(StrOutput{ka,2}));
    %count = fprintf(fid,'%20.6f',pos(1));
    %count = fprintf(fid,'%20.6f',pos(2));
    %count = fprintf(fid,'%20.6f\n',pos(3));
%end
%fclose(fid);


